package scu.maqiang.eigenvalue;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.Pair;
import scu.maqiang.numeric.SRMatrix;

public class SteklovEigenValue3DFine {

	public static void main(String[] args) {
		String dir = "SteklovEigenValue\\";
		Mesh3T4 mesh3D = new Mesh3T4().cube(64, 64, 64);
		//mesh3D.displayNodes();
		//mesh3D.displayElements();
		//mesh3D.displayBoundarys();
		int n = 8;
		double eps = 1.0/n;
		
		ScalarFunc chFun0 = (xy, label, param) -> xy[0]> 0.5 && xy[1] > 0.5 && xy[2] > 0.5? 1.0 : -1.0;
		mesh3D.removeElements(chFun0, null, 10);
		mesh3D.saveMesh_mesh(dir + "LMeshHomo128.mesh");
		mesh3D.toTecplot(dir+ "LMeshHomo128.dat");
		double[] param = {n};
		ScalarFunc chFun = (xy, label, para) -> {
			int nn = (int)para[0];
			double xCell = xy[0] * nn - (int)(xy[0] * nn);
			double yCell = xy[1] * nn - (int)(xy[1] * nn);
			double zCell = xy[2] * nn - (int)(xy[2] * nn);
			return (Math.abs(xCell- 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25 && Math.abs(zCell - 0.5) < 0.25)? 1.0:-1.0;	
		};
		mesh3D.removeElements(chFun, param, 8);
		mesh3D.toTecplot(dir+ "LMacroMesh64_8.dat");
		mesh3D.saveMesh_mesh(dir + "LMacroMesh3D64_8.mesh");
		
//		Mesh3T4 cellMesh3D = mg.cube(12, 12, 12);
//		param[0] = 1.0;
//		cellMesh3D.removeElements(chFun, param, 8);
//		cellMesh3D.saveMesh_mesh(dir + "CellMesh12.mesh");
		
		FES3T41 fs = new FES3T41(mesh3D);
		int dofHeat = fs.GetNdof();
		SRMatrix A = new SRMatrix(dofHeat);
		SRMatrix B = new SRMatrix(dofHeat);
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleHeatConvectSurface(new double[] {1.0}, BVPType.CONSIST_MASS, B, 8);
		double[] x = new double[dofHeat];
		Arrays.fill(x, 1.0);
		fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1, 2, 3, 4, 5, 6, 10);
		EigenSSolver solver = new EigenSSolver();
		
		int N = 10;
		double[] D = new double[N];
		double[][] V = new double[N][dofHeat];
		solver.inversePowerTran_PCGSSOR(A, B, x, 2.0, N, D, V);
		
		@SuppressWarnings("unchecked")
		Pair<Double, double[]>[] eigenpair = new Pair[N];
		for(int i = 0; i < N; i++) {
			eigenpair[i] = new Pair<Double, double[]>(D[i], V[i]);
		}
		
		Arrays.sort(eigenpair, (a, b) -> a.first > b.first ? 1 : -1);
		
		System.out.println(MVO.toString(Pair.getFirst(eigenpair)));
		double[] eiv = Pair.getFirst(eigenpair);
		String fileName = "eigenvalues128_8.txt";
		try {
			FileWriter fw = new FileWriter(fileName);
			for(int i = 0; i < N; i++) {
				fw.write((i + 1) + "\t" + eiv[i] + "\n");
			}
			fw.flush();
			fw.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
		
		mesh3D.toTecplot(dir + "SteklovEV3DL128_8.dat", Pair.getSecond(eigenpair));
	}

}
